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Abstract. 

We discuss results that have been obtained from the implementation of the 
initial round of testbeds for numerical relativity which was proposed in the first 
paper of the Apples with Apples Alliance. We present benchmark results for 
various codes which provide templates for analyzing the testbeds and to draw 
conclusions about various features of the codes. This allows us to sharpen the 
initial test specifications, design a new test and add theoretical insight. 



PACS numbers: 04.70.Bw, 04.25. Dm, 04.40.Nr, 98.80.Cq 



1. Introduction 

For decades, the field of numerical relativity has been dominated by an often painful 
quest for stable black-hole inspiral simulations. More than forty years after Hahn 
and Lindquist's first pioneering numerical simulation of colliding black holes J], this 
quest has recently turned into a gold-rush when Pretorius's breakthrough simulation 
[2] based on a harmonic code was followed by simultaneous invention of the "moving 
punctures" method by two independent groups [3l |4] . 

The primary motivation for solving the binary black hole problem in numerical 
relativity has however been to supply waveforms for gravitational wave detectors. 
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This goal demands an approach that goes beyond the efforts that have lead to an 
explosion in publications from the binary black hole community. Cross-validation of 
waveforms between different groups (and codes) and comparison with post-Newtonian 
predictions will be essential for numerical waveforms to be used in the computationally 
expensive searches conducted by the international gravitational wave community. The 
importance of cross-validation of numerical relativity results as a community effort was 
foreseen by the Apples with Apples Alliance (AwA) [5], which has presented a first 
round of standardized testbeds [6] . This first round comprises four tests with periodic 
boundaries, designed to efficiently exhibit code instability and inaccuracy. Instabilities 
currently receive less attention, since it has turned out that, paradoxically, binary 
black hole evolutions are in some sense a simpler problem than had been expected, 
and current codes evolving binary black holes do not typically show instabilities. 
The same codes will however have difficulties with some of the testbeds presented 
in the first round. The theoretical understanding of what works and what does not 
in numerical relativity is still very much an open problem. One crucial theoretical 
advance, which has been made since the publication of our first paper [6], is the 
development of a theory for well-posed second order in space, first order in time 
systems [U HI EJ [TOl HH [l2l US] , which has been extended to a basic understanding of 
numerical stability for such systems fTTJ [121 E] • 

Over the past years several groups have committed their test results to a 
publicly available data repository, with activities being coordinated via the web- 
site http://www.ApplesWithApples.org The purpose of the present paper is 



to document these developments and discuss their feedback with respect to code 
performance, to test improvement and to design further tests. While predating the 
binary black hole breakthroughs, we believe that the initial Apples with Apples tests 
and results are still valuable as providing a first testbed for a community effort in 
numerical relativity. 

The tests side-step many issues that would arise in a precise discussion of the 
binary black hole problem, such as the issue of boundaries. We make the natural 
choice of periodic boundaries for a first round of tests to isolate the performance 
of evolution algorithms. This is equivalent to evolution on the topology of a 3- 
torus in the absence of boundaries. However, in the context of general relativity, 
this introduces complications of a cosmological nature regarding the instability of 
Minkowski spacetime to perturbations on a compact manifold, as has been discussed 
in [6]. 

Establishing a paradigm for standardized testbeds for numerical relativity is a 
formidable task in itself. We can draw on experience from other fields, such as 
computational hydrodynamics where such testbeds have been used for a long time 
(for an overview of CFD testbed resources on the web, see e.g. [14]; for an example 
of initial value ordinary differential equation (ODE) test-suites see [15]). However, 
general relativity comes with its own issues that introduce extra complications. First 
of all, it is important to realize that the numerical relativity community is small, 
with very limited available manpower. In contrast to the size of the field, we are 
trying to solve many difficult problems at the same time. Numerical methods are 
being developed in parallel with the formulation of the continuum problem, with the 
construction of physically relevant initial data sets and with the unraveling of the 
physical processes involved in the systems under investigation. All of this is, so far, 
without the help of comparison with experiments. Groups working in the field are 
faced with many fundamental questions in designing their approaches. Codes are in 
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a state of flux that makes careful documentation easy to postpone. A good example 
is the issue of boundaries, which can be taken to be either a cubic grid boundary or 
a smooth spherical boundary, which can either be mapped to infinity or given some 
finite artificial location, and which are further complicated by gauge freedom and 
the requirements of constraint preservation. Useful comparison of the wide variety of 
resulting codes requires simple tests which isolate an important facet of the problem. 

We distinguish two fundamentally different types of testbed: The first type 
compares different codes and methods in the treatment of a physically interesting set 
of solutions. In the context of the binary black hole problem, a detailed comparison 
of nonspinning equal-mass inspiral would be a natural example. The second type are 
idealized situations, such as the "shock tube test" [l6j in computational fluid dynamics. 
This is the type of testbed we discuss in the present paper, where we restrict ourselves 
to a greatly simplified first set of tests [6]: periodic grids and strict test specifications, 
which as far as practicable define all the details of a simulation except the formulation 
of the Einstein equations. Our experience with the first round of testbeds confirms 
this decision: even the analysis of these simple situations has proved quite challenging. 
Our conclusions in Sec. [8] discuss how the experience from the present round of tests 
can be used in our development of black hole tests. 

We identify five main aims of standardized tests of the "idealized" type: 

(i) Standardized tests should provide the young and fast-changing community of 
numerical relativists with a common reference frame which will help integrate 
different efforts to produce a coherent picture of what works and what does not, 
and thus reduce the dependence on anecdote and fashion. 

(ii) Tests should be efficient in revealing instabilities or other weaknesses of an 
algorithm, both regarding simplicity of the analysis, run time and implementation. 

(iii) Tests should help identify where problems come from, as a step toward 
improvement of the algorithms. 

(iv) Tests should facilitate comparisons between approaches regarding different 
continuum formulations, spatial discretizations, time integrators, uses of artificial 
dissipation, etc. 

(v) The development of testbeds should eventually lead to useful code comparisons 
for judging the validity of physically interesting simulations, e.g. the binary black 
hole problem. 

Point (i) has been addressed by organizing this project as a community initiative, 
which seeks broad participation and provides test results via web pages and a 
CVS repository [5]. Regarding point (ii), in this paper we review our original test 
specifications and propose modifications to promote efficiency. Point (iii) is essential 
for the character of this paper: we focus on presenting test results as a template for 
analyzing and interpreting results, rather than just presenting the broadest possible 
listing of test output for a maximal number of codes. We feel that it is essential to 
stress this point: tests which do not directly correspond to a physically interesting 
situation are only valuable if they improve our understanding of what really goes 
on with a certain code. Only then can we hope to carry over test benefits to other 
situations. Such analysis does of course require a certain effort. 

Point (iv) is dealt with by providing "standard candle results" in the CVS 
repository, i.e., benchmarks that have been obtained with very strictly defined 
specifications. Point (v) represents the ultimate goal of the AwA Alliance. 
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The analysis of test results has led to better understanding of the four original 
standardized tests and has led to some improvements in their specifications. We also 
have added a new shifted gauge wave test, which closes a gap regarding the ability of 
a code to handle a shift. The revised specifications for the five tests are detailed in 
[Appendix A[ The major changes from the specifications in [6] are 

• Robust stability test: The rules for how the data should scale with resolution 
have been changed; the criteria for passing the test, has been restated. 

• Linearized wave test: No changes. 

• Gauge wave test: The original tests amplitudes A = .01 and A — .1 have been 
replaced with A = .5. 

• Shifted gauge wave test: This new test has been added. 

• Gowdy wave test: No changes. 

We have also dropped the original requiremnt that the tests be run with a iterative 
Crank-Nicholson integrator. Conclusions from the test results and our experiences 
with the testing procedures, along with the reasons behind the changes and additions 
in the standard tests, are summarized in Sec. [8] 

The code descriptions and test data on which this paper is based are described 
in Sec. |2l The results for the original four standardized tests are discussed in Sees. |3l 
m \5\ and |7l Discussion of the shifted gauge wave test and some benchmarks are given 
in Sec. HI 

The plots presented in this paper are based upon test output in the CVS 
repository. Many of these tests were run with codes in which artificial dissipation 
was only introduced implicitly through the use an iterated Crank-Nicholson (ICN) 
time integrator. It had been a naive hope at the beginning of this project that the 
use of ICN might provide a way to standardize the introduction of dissipation. Most 
numerical relativity groups now use Runge-Kutta time integrators with the explicit 
addition of Kreiss-Oliger dissipation (see [Appendix C."2 ). It has been found that many 



of the test results presented here could be greatly improved by such explicit use of 
dissipation. In addition to artificial dissipation, most codes that simulate binary black 
holes use higher order approximations than the second order accurate codes being 
compared here. Consequently, we want to emphasize that the results exhibited in this 
paper should not be used to make judgments on particular approaches, but that our 
purpose is to assess and improve the test suite and to provide a basis for future code 
comparisons. 



2. Code descriptions 

In order to ensure a consistent presentation of test output, we present a brief account 
of the numerical codes and algorithms which have been used to produce the data on 
which this paper is based. All data are publicly available via the CVS repository 
(see [S| for details). The four original standardized tests are denoted by ROBUST 
(the Robust Stabihty Test), LINEAR (the Linear Wave Test), GAUGE (the Gauge 
Wave Test) and GOWDY (the Gowdy Wave Test). Table 1 summarizes the output 
data that have been submitted for the various codes. 

The usefulness of this data depends upon good code documentation. It is beyond 
the scope of this paper to provide such documentation for all the codes involved. 
However, we will outline some basic code information which is necessary to interpret 
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CODE 


ROBUST 


LINEAR 


GAUGE 


GOWDY 


AbigcLharm 


++ 


++ 


++ 


++ 


AELCactusEinstcinADM 


+ 






++ 


Kranc_FrccADM 


+ 


+ 


+ 


+ 


CCATIE.BSSN 


++ 


++ 


++ 


++ 


Kranc3SSN 


++ 


++ 


++ 


++ 


LazEv_BSSN 


++ 


++ 


++ 


++ 


HarmNaivc 


++ 


++ 


++ 


++ 


KrancNOR 


++ 


++ 


++ 


++ 


KrancFN 


++ 


++ 


++ 




LSU.HypcrGR 


++ 


++ 


++ 


++ 



Table 1. Test output and codes considered in this article. The code abbreviations 
are explained below, along with a description of the finite difference algorithm. A 
"++" indicates a full complement of test output in the CVS, a "+ indicates partial 
output which has been used for our analysis, a "— " indicates partial output on 
which no meaningful conclusions could be drawn and a " " indicates no output. 



the test results. The complexity of this task is somewhat alleviated because all the 
codes represented here follow a method of lines approach. We will organize the code 
descriptions along the following guidelines. 

• A description of the continuum formulation^ including a list of all variables, their 
associated evolution equations and constraints (both differential and algebraic), 
equations governing the lapse and shift and a specification of any free parameters. 
Terms and differential operators in the equations should be ordered in the way 
that they are approximated by finite difference expressions in order to avoid 
ambiguities associated with the Leibniz rule. The hyperbolicity classification 
should be provided, if known. 

• A description of the semi-discrete system, describing the spatial finite difference 
equations on each time level, including the rules for discretizing partial derivatives 
as centered or one-sided finite differences and any other discretization techniques, 
such as spatial averaging or dissipation. For complicated systems, the finite 
difference rules may be specified only for the principal part, with further 
details supplied by references. (Here we provide some basic reference material 
in Appendix B| and [Appendix C| for compactness of presentation.) 



• A description of the numerical time update scheme. All manipulations of 
data between intermediate time steps should be specified, such as enforcing a 
constraint. 

As an example, we consider two inequivalent algorithms for the wave equation 
D(j) = (with unit lapse, zero shift and spatial metric 7ij), which should be expected 
to result in different code performance. In both cases the second order in time system 
is reduced to first order in time by introducing the variable n = dt4>, and applying, 
say, 4th order Runge-Kutta (see |Appendix C[ ) to the ODEs of the semi-discrete system 
obtained using the method of lines. Two different codes can based upon the following 
descriptions. 

Description I: 
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(i) The continuum system is 

9*0 = TT, (1) 

dtn = ±d,{^f^d,<t>). (2) 

(ii) The semi-discrete version is obtained by replacing all partial derivatives in ([2]) by 
centered differences: 

where Doi is the centered difference operator Dg applied in direction i 
(see [Appendix C.l[ ) . 

Description II (incquivalent with I): 

(i) The continuum system is 

9*0 = TT, (3) 

dtiT = 7''a«9,0 + J-d,{^Y^)d,cf>. (4) 

(ii) The semi-discrete version is obtained by replacing the partial derivatives in ([2]) 
by centered differences according to 

= YW+,D^,cf> + ^D^,{^f^)D^,<ly (5) 

where -D+i and D^i represent forward and backward centered finite differences 
in the respective directions (see [Appendix C.T ). 

The codes resulting from these two descriptions produce substantially different 
performance because of the "checkerboard" design of the stencil used in description 
I. Descriptions of the specific codes used in this paper are given in [Appendix B[ 



3. Robust stability test 



The robust stability test was intended as a first screen to eliminate many unstable 
evolution algorithms. The particular importance of this test was due to the fact that 
instabilities of numerical codes appeared as a prime obstacle to "solve" the binary black 
hole problem, and essentially no theoretical understanding was available to discuss 
the well-posedness and numerical stability of first order in time, second order in space 
formulations of the Einstein equations, which have been and still are popular in the 
field. Recently, a theoretical framework has become available to discuss the well- 
posedness and numerical stability of such mixed order formulations of the Einstein 
equations [17llll|8l|9l|10l [HI dH (El [13] , and it has been extended to the problem 
of discretizing the equations in the context of the method of lines [TTl [T^l [13] • As a 
consequence of both the recent breakthroughs in the binary black hole problem and 
the theoretical advances, numerical stability has become a relatively minor issue in 
practice (although there certainly remain interesting mathematical questions to be 
pursued). We thus restrict ourselves to a minimal discussion here, as is sufficient 
to understand the data available in our test results repository. For a more in-depth 
discussion of theoretical and practical aspects of numerical stability and the robust 



Implementation of standard testbeds for numerical relativity 



7 



stability test we refer to , which has been directly motivated by numerical results 
obtained within this project. 

While the other tests give quantitative information about an evolution system, 
e.g. the magnitude of the numerical error, the result of the robust stability test is 
"pass" or "fail" . A stable numerical algorithm is only possible if the underlying 
continuum problem is well-posed |19j . In the well-posed case an instability might 
still arise, either from the numerical technique or from the existence of an exponential 
mode in the continuum problem. The test is designed to avoid continuum instabilities 
by considering small perturbations of the Minkowski metric. In addition to providing 
efficient detection of unstable numerical algorithms (or coding errors) affecting the 
principal part of the evolution system, it is also intended to spot instabilities arising 
from ill-posed systems, such as weakly hyperbolic systems. 

As an example, consider the weakly hyperbolic system 

v,t = (6) 
with the periodic solutions 

u = Lot cos Lj{t x) ^ f = sin (t 4- x) (7) 

Lo = 27rTO , m ~ 1, 2, 3, ... 
on the domain —.5 < x < .5. In terms of the L2 norm 

iV = y {u^ + v^)dx^ ^ , (8) 

the Cauchy data for d?]) at i = 0, 

u = 0, w = sinwx, (9) 

has norm iV(0) = l/\/2. However, because of N{t) ^ cut for large uj. This leads 
to a violation of the well-posedness requirement that in any finite time interval 

N{t) < yle^*Ar(0), (10) 

in terms of constants A and K independent of the Cauchy data. 

For discretized systems we can not test well-posedness directly, but rather we test 
the analogous concept of numerical stability, i.e., we aim at establishing the existence 
of constants A and K, which give rise to the bound 



where is the solution of the discrete system at time t„ ~ nk. The test is passed if 
such a bound can be established, and is failed otherwise. In the discretized version of a 
weakly hyperbolic problem, with grid displacement h, the perturbation of a simulation 
by random initial data can be expected to excite numerical error which grows linearly 
in time according to u ~ t/h, corresponding to the shortest wave number uj l/h. 
This would then lead to secular error growth which increases with resolution. Although 
the system (|6]) is well-posed with respect to a stronger norm including a term, 
a generic perturbation of ^ by lower order terms would nevertheless produce an 
exponentially growing instability which cannot be bounded. See |20j for a more general 
discussion of such weakly hyperbolic systems. 

The key idea of setting initial data for this test is to distribute energy roughly 
equally over all frequencies. This is a particularly efficient way to reveal growing 
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modes if the growth rate increases with resolution, as is the case if the discretization 
is unstable or if the continuum problem is ill-posed. In our test we use a spectrum 
generated by random initial data. 

The robust stability test as formulated here tests numerical stability in the linear, 
constant coefficient regime. It is based upon small random perturbations of Minkowski 
space, with the initial data consisting of random numbers e applied at each grid point 
to every code variable requiring initialization. In numerical evolution, where machine 
precision takes the place of e, a code that cannot stably evolve such random noise 
would be unable to evolve smooth initial data. 

In spite of its simplicity, our experience has shown that the robust stability 
test exhibits various subtle difficulties in designing a single test prescription that is 
universally effective for all evolution systems and numerical methods. Some particular 
problems are: 

• For random initial data, where a significant part of the total energy is in 
high frequencies, dissipation has a large effect. Some intrinsic dissipation 
is unavoidable in finite difference evolution algorithms, and adding artificial 
dissipation may be necessary to stabilize certain algorithms jllj . and insufficient 
to stabilize others (such as algorithms for weakly hyperbolic systems). 
Simulations of variable coefficient, nonlinear systems normally require numerical 
dissipation to obtain a stable evolution, e.g. by adding Kreiss-Oliger type 
dissipation [20] (see [Appendix C3 ). Dissipation can however increase the time 



scale on which instabilities become apparent. The detailed way dissipation affects 
instabilities varies with the spatial discretization (we only consider second order 
approximations here), with the time integrator, with the grid resolution and with 
the Courant number. 

• As discussed in the above example, well-posedness and numerical stability are 
defined with respect to a certain norm. Using an inappropriate norm can yield 
misleading results. Second order systems require different norms than first order 
systems pT] . 

• Numerical stability of an explicit time integration algorithm can only be expected 
if the time step is appropriately restricted by a Courant-Friedrichs-Lewy (CFL) 
condition. It is important to distinguish between resolution dependent blowup 
associated with ill-posedness from blowup resulting from a CFL violation. For 
sufficiently complicated 3D algorithms, the CFL limit might not be readily 
deduced from analytic arguments. As an example, exponential growth of the 
ADM algorithm was mistakenly provided as an illustration of a failed robust 
stability test in [6] . It took subsequent testing and analysis to reveal that this 
exponential growth resulted from a CFL violation and that otherwise the weakly 
hyperbolic instability of ADM resulted in a secular (linear in time) growth. 

As a result of such considerations, we will not try to present a single universally 
applicable specification for the robust stability test. Instead, while keeping the original 
spirit of the test as a simple and useful first screen, we propose some changes in the 
guidelines, as discussed below. 

An important issue when performing stability tests is whether the high frequency 
modes are damped. This has important bearing on the long-time behavior of the 
robust stability test: all damped modes will decay in time; eventually the undamped 
frequencies of the discrete system will dominate the signal. If an analysis of damping 
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factors has not been performed, the test can therefore also be useful in detecting 
the spectrum of frequencies which are not damped. It has been pointed out in [TT] 
for standard discretizations of first order in space systems that the "checkerboard" 
mode is undamped, while for typical second order systems it is damped. Since the 
"checkerboard" mode is not realized on grids with an odd number of points, we adopt 
the practice of always using an even number of grid points so as not to muzzle such a 
potential instability. 

In our original specifications, we proposed the relatively large time step dt = 
O.bdx, which turned out to be larger than the CFL limit for the ADM system. Since 
a smaller dt also decreases the amount of dissipation inherent in a time integrator, 
we now propose a relatively small time step to avoid distortion of results due to 
dissipation. Common time integrators in current practice in numerical relativity 
are ICN, RK3 and RK4 (sorted by decreasing internal amount of dissipation). A 
sufficiently small time step would yield similar results for all of them. We therefore 
propose to run with dt = Q.ldx, which can be further reduced in case of doubt. See 
[Appendix A.l for details. 

For systems that use variables which correspond to spatial derivatives of the ADM 
3-metric and extrinsic curvature, an ambiguity arises: noise can be added uniformly to 
all variables, or to the ADM initial data before taking derivatives. There are similar 
ambiguities in second order systems regarding how the range of the random numbers 
should scale with resolution. For uniformity of description, we propose to do the 
simplest thing, namely to apply noise to all evolution variables in the same way. We 
propose the range of ±10""'^° for all variables, the same range used for the lowest 
resolution in the original specifications. 

Following common practice at the time, the Hamiltonian constraint was used to 
analyze test results. Again following [llj . we now propose a pass/fail analysis based 
upon whether the time behavior of the norm satisfies (jlip. 

Our core test specification combines both ID and 3D features by running in a 
thin channel along the x-axis. The use of 4 distinct gridpoints in the y and z directions 
allows for the checkerboard mode (ghost points may be necessary depending upon the 
numerical scheme). The generalization to a full cube 3D test is straightforward, and 
may add further clarification in case of dubious results. 

The test should be run until one is confident that dissipation effects do not cloud 
the result. Without artificial dissipation, a runtime of one crossing time, using output 
at every time step, is usually sufficient. This corresponds to 500p time steps, for 
a given resolution p (see [Appendix A} . The test is passed if the norm satisfies the 
inequality pT|) for all resolutions, for a fixed choice of A and K. 

Instabilities caused by the ill-posedness of the evolution system (or by coding 
errors in treating the principal part), are already apparent in one-dimensional tests, 
which can be performed quickly and economically. An example of how this analysis 
works is given in Fig. [TJ The way that the slope of the error vs time depends 
upon resolution shows that the AbigeLharm code, which is based upon a symmetric 
hyperbolic formulation, passes the test; whereas the HarmNaive code, which is based 
upon a weakly hyperbolic formulation, fails the test. 



4. Linearized wave test 



A prime physical objective of numerical relativity is to compute the waveform from 
a system of black holes and neutron stars. This test checks the ability of a code to 



Implementation of standard testbeds for numerical relativity 



10 




Figure 1. Convergence results for the robust stability test with the AbigeLharm 
(left) and HarmNaive (right) codes, for runs of 1 crossing time. The graphs show 
the error in g^x as a function of time, obtained by subtracting 1 from its L2 norm. 
As seen from the slopes of the graphs, the Abigeljiarm code (left) passes the test, 
because there is no increasing rate of error growth with higher resolution p, while 
the HarmNaive code (right) fails the test because the growth rate increases with 
resolution. 



propagate a linearized gravitational wave, which is a minimally necessary attribute 
for reliable wave extraction from strong sources. Test specifications are given in 
[Appendix A.2| 

The test checks the accuracy of the code in propagating both the amplitude and 
phase of the wave. It can reveal whether excessive dissipation has been necessary for 
good long term performance in the robust stability test. For the p = 1 coarsest grid 
(TV = 50 grid zones), there is not enough resolution for second order accurate codes to 
obtain accurate phase propagation and the corresponding runs should only be viewed 
as an economical first check on the code. The most useful comparisons are with the 
p = 4 grid. 

Fig. [2] compares snapshots of the ID wave after 1000 crossing times which were 
obtained with a variety of codes using the p = 4 finest grid. For reference, the 
exact waveform is also plotted. The snapshots for three of the codes, AbigeLharm, 
HarmNaive and LazEv_BSSN, are very similar and provide a good benchmark for the 
accuracy that can be achieved at this resolution. They very closely match the exact 
solution in amplitude but show a phase delay, similar to the delay seen in the following 
gauge wave test. It should be expected that phase accuracy could be improved by 
going to fourth order accurate methods. Some snapshots of the corresponding error 
are displayed in Fig. [3] Except for the two codes with the largest phase error, the 
error at 1000 crossing times is confined to a small band. By monitoring the growth 
of the error during the evolution, it was verified that no overall multiple of 27r phase 
error is concealed in the snapshots of Fig. [21 

In addition, the plots of the Hamiltonian in Fig. [J show no rapidly growing 
constraint violating instabilities in this linear regime. The secular instability of 
Harm_Naive, which was discussed in the robust stability test, is evident but it does 
not introduce a large error in this test. This illustrates that instabilities associated 
with a weakly hyperbolic system arc not necessarily evident in linearized tests where, 
as discussed in Sec. [3j the unstable modes only grow secularly in time. The KrancFN 
code gives good accuracy for the amplitude but a much larger error in phase. The 
CCATIE code shows poor accuracy in both phase and error. It is beyond the scope 
of this paper to explain the discrepancy between the performance of the two BSSN 
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Figure 2. Comparison snapsliots of gyy{x) — 1 at t = 1000 for tlie ID linearized 
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Figure 3. Comparison snapshots of the error £ in gyy{x) at t = 1000 for the ID 
linearized wave test, with p = 4 resolution. 



codes. 

The ID linear wave test is simple and economical to perform. Although the test 
is not very demanding, the results for the metric component gyy in Figs. [2] and [3] 
show that it provides a benchmark which can be useful to identify weaknesses in code 
performance. The 2D tests require more computer time and the results were typically 
in line with expectations from the ID results. 



5. Gauge v^rave test 

The gauge wave test is based on a nonlinear gauge transformation of Minkowski 
spacetime. Although the correct solution is a flat spacetime, nonlinear effects and 
the nontrivial geometry of the time slices can easily trigger continuum instabilities 
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Figure 4. Comparison of time dependence of the L^a norm of the Hamiltonian 
constraint WHW, shown on a logarithmic scale, for the ID linearized wave test with 
p = 4 resolution. 



in the equations. For simple examples of such effects see |21| for a nonlinear -wave 
equation on flat space, designed to model problems arising in this testbed, and |22| 
for a linear example of ho-w nontrivial geometry of the slicing can trigger instabilities 
already for the Max-well equations. 

Our original specifications [6] "were to run the test "with amplitudes A = 0.01 and 
A = 0.1. Many codes have been sufRciently improved to handle larger amplitudes, 
"which is generally more efficient in detecting instabilities "with smaller run times. 
Accordingly, -we specify an amplitude of ^ = 0.5 in the revised test details given 
in [Appendix A. 3 

While the gauge "wave metric has a rather simple form, the test proved to be 
challenging for most evolution codes. One anticipated source of gro-wing error is the 
instability of a flat space -with topology [6j. Another problem is the existence 
of a family of harmonic, exponential gauge modes corresponding to the substitution 
H e^^H (for arbitrary A) in the metric (|A.10[) [H]. The testbed itself corresponds to 
A = 0, but numerical error can easily excite this mode and lead to exponential gro-wth 
of the -wave amplitude. Other instabilities may be present in individual systems, 
depending on the detailed form of the reduced evolution system for the particular 
formulation. Some of these instabilities can be identified by looking at the gro-wth 
of the constraints for the formulation. In addition to instabilities that correspond 
to solutions of the continuum problem, individual codes may suffer from numerical 
instabilities depending on the discretization schemes. These "would typically be seen 
as high frequency modes and, for -well-posed systems, can be cured by adding artificial 
dissipation to the numerical algorithm. 

Figure [S] sho"ws the time evolution of the Hamiltonian constraint for the various 
codes. The negligible violation of the Hamiltonian constraint by the harmonic codes 
can be attributed to the fact that the harmonic coordinate conditions are used to 
shift the role of the constraint to an evolution equation. Note that the BSSN codes 
sho"w rapid gro-wth of Hamiltonian constraint violation. So far no BSSN code has 
demonstrated satisfactory performance for this test, and for brevity -we do not include 
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Figure 5. Time dependence of the Loo norm of the Hamiltonian constraint WHW, 
shown on a logarithmic scale, for the ID gauge wave test with resolution p = 4 
and amplitude A = 0.1. 



5.1. Results 

5.1.1. Results for the AbigeLharm Code For this particular testbed most components 
of the densitized metric g'^'^ ~ \/~99^'^ have trivial values, the non-trivial ones being 

-gyy = 5- = H (12) 

The original implementation of the Abigel code based upon (jB.9p leads to a numerically 
stable and convergent code, with no high frequency modes generated. However, as 
shown by the dramatic growth of the rescaled error plotted in Fig. [51 the gauge wave 
excites exponential modes g^^ — ~ 1 — e'*'*iJ, A > 0. This can be understood [2T|in 
terms of solutions of the harmonic system whose densitized metric components arc all 
trivial except for 

r' = r = F{t,x). (13) 
The resulting source term S*^" in (jB.9p vanishes except for the components 



gyy — 



F 



(14) 



The PDE for F{t,x), which results from inserting (fT3|) into (jB.9p . reduces to 
{—df + d^)\ogF = 0, which admits the exponential solutions F = e^*H. These 

(15:91). so 



solutions satisfy the harmonic constraints and the reduced harmonic system 
that they are also solutions of the full Einstein equations. Therefore all codes using 
harmonic gauge conditions might be expected to excite this mode. 

In the case of the AbigeLharm code, these modes were suppressed by building 
semi-discrete conservation laws into the code which, for the gauge wave initial data, 
would not be obeyed by the exponential solution. Namely, by writing (|B.9[) in the flux- 
conservative form (|B.10p . the principle part of the resulting equation has vanishing 
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Figure 6. Comparison of code performance between the non-fiux-conservative 
(non-FC) and flux-conservative (FC) versions of the AbigcLharm code, showing 
graphs of g~"{x) at t = 100 for a gauge wave of ampHtude A = 0.5 on the 
p = 2 grid. In the non-FC case the graph is rescaled by the average of the 
plotted function, showing g^^/avg{g") g^^/ exp(29.8). The good overlap of 
this rescaled function with the analytic value clearly indicates that the dominant 
error of the non-FC code is a multiplicative function of t. Measurements at 
t = 100 for the non-FC code show that logarithm of the spatial average of 
scales roughly as {dx)'^ , i.e., log(ai)g(g^^)p=i) 110.8, log(av(;(g^^ )p=:2) ~ 
29.8, log(ai'(;(g'''')p=4) as 7.52, suggesting that the multiplicative error has 
exponential growth of the form cxp(0(((ia:)^) ■ t). 



source term, S^'^ = 0, for this test. A summation by parts numerical algorithm then 
gives rise to the semi-discrete conservation law 



While this is a non-generic result (most space-times would give a non-zero source 
term), building this conservation law into the principal part of the system has proved 
effective not only in this particular case but in the other Apples with Apples tests 
considered in this paper, as well as in further proposed tests [HI [531 [H]. 

As shown in Figs. [71and[51 the flux-conservative code does not develop exponential 
error modes when running with the original ICN integrator (see [23] for similar results 
with RK4.) The main source of error is phase error which converges to zero as the 
grid is refined. In order to further illustrate this point. Figs. [3 and [3 give test results 
for both the ID and 2D versions with amplitudes of A = 0.01, 0.1, 0.5. 

5.1.2. Results for the HarmNaive System This naive harmonic system, although 
weakly hyperbolic, behaves identical to the symmetric hyperbolic AbigclJiarm code 
for this testbed. This can be understood given that the RHS for the mixed space-time 
components of the evolution system vanish, i.e. 



which implies that the time-time component of the RHS also vanishes, i.e., 




(15) 



I,J,K 



dtg'* = -djf^ = 0, 



(16) 



dtg 



■tt 



djf' = 0. 



(17) 
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Figure 7. Convergence results for the ID gauge wave simulation writh the 
Abigeljiarm code, for amplitudes of A = 0.01 (left) and A = 0.1 (right). The 
graphs show the Loo norm of the error in Qxx, defined as g^^^ = g^x"^ ~ 9xx'^ ^ 
a function of time, and rescaled by a factor of 1/p^ . As seen from the graphs, the 
lower amplitude runs give no new information. 
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Figure 8. Convergence results for the 2D gauge wave simulation with the 
Abigeljiarm code, for amplitude A = 0.5. The left graph shows the Lao norm of 
the error in gxx, rescaled by a factor of 1/p^, as a function of time; while the right 
graph shows the same rescaled error norm for the violation of the Hamiltonian 
constraint H. For the AbigeLharm code, the vanishing of the Hamiltonian 
constraint is an algebraic identity, making H of order roundoff. As a result, 
the constraint violation is super convergent. The lower amplitude runs revealed 
no new features. 

The test-results confirm this. 

As expected, tests for the ADM-system also behave identically, since the naive 
harmonic system can be understood as a formulation of the ADM-system in the 
harmonic gauge. We therefore skip a separate discussion of the ADM-system. 

5. 2. Results for the KrancFN and KrancNOR Systems 

Besides the harmonic codes, KrancFN was the only other code that was able to run 
for 1000 crossing times for an amplitude A= .1. At the end of the run, Fig. [S] shows 
that long wavelength growth due to the e^*i/ instability of the wave amplitude has 
become appreciable. 

The KrancNOR code picks up the e^*H instability at a faster rate and, although 
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it shows clear 2nd order convergent at early times, it crashes at < w 44. The snapshot 
in Fig. [5] shows that the error at the end of the run is almost exactly in the e^^H 
mode. 
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Figure 9. Comparison snapshots of for tlie ID gauge wave witii amplitude 

A = 0.1 at the end of a run with p = 4 resolution. For the Abigel and KrancFN 
codes, the run lasts the full 1000 crossing times. The KrancNOR code crashes at 
t = 44. 



6. Shifted gauge wave testbed 

In the shifted gauge wave test (jA.lSp we have identified two types of instability 
One, which is analogous to the instability of the gauge wave, arises from the A- 
parameter family of vacuum metrics 

dsl = e^\-dt^ + dx^) + dy"^ + dz^ + Hkakfjdx^dx^ , (18) 

which reduces to the shifted gauge wave for A = 0. The other is an instability peculiar 
to harmonic (or generalized harmonic) evolution codes, where the Einstein equations 
are satisfied only indirectly through the harmonic conditions. The metric 

dsl = -dt^ + dx^ + dy"^ + dz^ + [h - I + e/*^ ko.kpdx'^dx^ , (19) 

where 

Ad (2TT{x-t) 



i = t--cos^^^j, (20) 

satisfies the reduced harmonic evolution equations (|B.9p . The simulation of the shifted 
gauge wave by any evolution code based upon a standard reduction of Einstein's 
equations to harmonic form can be expected to excite this instability. 

The test was developed in conjunction with the AbigeLharm code [23j . For ID 
runs with the /C = 4 resolution, it was found that the evolution equation (jB.lOp excited 
the instability (|19|) on a timescalc t w 500. Further investigation showed that this 
instability could be suppressed by adjusting (jB.10[) according to 

^ g,.u _ j^^u^ ^21) 
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where A^" = when the harmonic constraints 

:= ^(S^g'"' -H'') = (22) 

V-9 

are satisfied. Particularly effective were the constraint adjustments 

^""-^^^C^C^ fo>0, (23) 

where Cpcr is the natural metric of signature (+ + ++) associated with the Cauchy 
slicing, and 

A^' = --^C"do.{^gg''''), c>0. (24) 

This is exhibited in Fig. (TUl which shows for a run with amplitude A = 0.5 that these 
constraint adjustments suppress instabilities for the entire 1000 crossing time duration 
of the test. 




Figure 10. Plots of the Loo error E(t) in gxx obtained with the Abigel code 
for the ID shifted gauge wave test with ampUtude A = .5 and resolution p = 4. 
Results are compared for the constraint adjustment II23I I with b = 1, the constraint 
adjustment I I24I I with c = 1 and the bare algorithm. The two adjustments show 
very similar error and both give excellent suppression of the unstable mode excited 
by the bare algorithm. 

Results for the shifted gauge wave tests are also available from the CVS repository 
for BSSN codes. In this case, as in the standard gauge wave test, results are not 
satisfactory, and suggest further analysis, which is beyond the scope of this paper. 
Results obtained with the Kranc_BSSN code and a very small value of the dissipation 
parameter (a = 0.001, see Eq. ()C.7| ) for the medium amplitude A — 0.1 are shown 
in Fig. 111! While the code shows second order convergence for several crossing times, 
rather quickly an instability develops that eventually crashes the code. As expected, 
the instability develops slower for the lower amplitude A = 0.01, and faster for A = 0.5, 
where the code crashes within roughly one crossing time. Similar results are also 
available in the CVS repository for the CCATIE code. 

Results for the shifted gauge wave test have also been obtained [3S] using the 
Caltech-Corncll group's spectral version of a code based upon the Kidder-Scheel- 
Teukolsky formulation of the Einstein equations [26]. For the ID test with A = .5, 
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Figure 11. Performance of Kranc_BSSN for the shifted gauge wave with 
amphtude A = 0.1 and a dissipation value of ct = 0.001. Left paneh The L2-norm 
of the Hamiltonian constraint plotted vs. time for resolutions p = 1,2,4 (short- 
dashed, long-dashed, full line). Right panel: Convergence test for the L2-norm 
of trK for resolutions p = 1,2,4. Second order convergence is lost after a few 
crossing times. 



they encountered nonlinear instabilities associated with aliasing after a few crossing 
times. There are standard filtering techniques to deal with such aliasing problems. 
By filtering the top 1/3 spectral coefficients, they found that the evolutions could be 
extended as far as i = 60, but further improvements by filtering did not seem possible. 
The group has not yet reported results for their current spectral code which is based 
upon a generalized harmonic formulation. 

7. Gowdy wave test 

The previous tests involve spacetimcs with small curvature. The Gowdy wave test 
is based upon a strongly curved exact solution for an expanding vacuum universe 
containing a plane polarized gravitational wave propagating around a 3-torus |27j . 
See |28| for a recent review. The metric has the form 

ds^ = i-i/2eV2(_rfi2 ^ ^^2) ^ ^(gP^^2 ^ e-Prfy2)^ (25) 

where P(t, z) and A(t, z) depend periodically on z and the time coordinate t increases 
as the universe expands, with a cosmological type singularity at i = 0. The detailed 
tests specifications given in [Appendix A. 5] were designed so that neither very large nor 
very small numbers enter in the initial data. 

In the expanding direction, the qualitative behavior of the solution is 
characterized by P slowly decaying to zero while A grows linearly, with both P and A 
exhibiting gravitational wave oscillations. The linear growth of A leads to exponential 
growth of gzz, so that code accuracy is tested in a harsh situation. This makes 
evolution with a 3D code difficult compared with the direct ID evolution of P used 
in numerical studies of the approach to the cosmological singularity [52] 

The performance of the various codes in the expanding direction is illustrated by 
the output for the trace of the extrinsic curvature K shown in Fig. [121 Although not 
apparent from the figure, the HarmNaive code crashes abruptly at i = 8, as might 
be expected of a weakly hyperbolic system in the nonlinear regime. Even though 
the analytic value of K is negative and asymptotes to zero with the expansion, short 
wavelength error in the LS_HyperGR and LazEv_BSSN codes triggers an instability 
leading to a collapsing mode with K > 0. This is illustrated for the LazEvJBSSN run 
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in the snapshot of Fig.[T31 which shows the error in gzz{t, x) a.t t = 13 just before the 
run crashes. The superposition of short wavelength error with the long wavelength 
truncation error from the signal is evident. 
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Figure 12. Comparison plots of the trace of the extrinsic curvature K for 
the polarized Gowdy wave evolved in the expanding direction with the p = 4 
resolution. Analytically K is spatially homogeneous; the plots show its maximum 
value over the numerical grid. 
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Figure 13. Plot of the error £{z) in gzz for the polarized Gowdy Wave evolved 
in the expanding direction with p = 4 resolution with the 2"'*order accurate 
LazEv_BSSN code. The error, plotted at t = 13 just before the code crashes, 
shows a large short wavelength component which can be controlled by dissipation. 

Further experiments with the LazEvJBSSN code showed that this short 
wavelength instability could be controlled by numerical dissipation and that 
the accuracy could be further improved by using fourth order finite difference 
approximations. For the expanding Gowdy test, this is illustrated in the plots of the 
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rescaled error in the left portion of Fig. [TJ] which indicate fourth order convergence. 
However, the error stiU exhibits poor long term accuracy. In the right portion of 
Fig. I14[ we also display the error in the second order accurate AbigeLHarm code. 
Both the second order and fourth order codes have approximately the same long term 
rate of growth due to the underlying exponential growth in gzz- 




Figure 14. Convergence results for the norm of the error [[^(i)!! in g^^ 
(logarithmic scale) for the polarized Gowdy wave evolved in the expanding 
direction. On the left, the results for the p = 2 resolution have been rescaled 
by 1/16 for the 4"* order accurate LazEv_BSSN code with dissipation. The 
results indicate stability and convergence but do not give long term accuracy. 
On the right, the error for the p = 2 resolution has been rescaled by 1/4 for 
the 2"'' order system Abigel_Harm code, again showing stability and convergence. 
Both codes exhibit roughly the same long term rate of error growth expected from 
the exponential growth of g^z- 

The Gowdy test is run in both future and past time directions because analytical 
studies [30] and numerical experiments [5^ [31] indicate that the sign of the extrinsic 
curvature may have important consequences for constraint violation. The subsidiary 
system governing constraint propagation can lead to unstable departure from the 
constraint hypcrsurface. As an example, in a hyperboloidal slicing of Minkowski space 
with unit lapse and zero shift, the electromagnetic constraint C = VaE"^ satisfies 
C(t) = C(0)e^* when the standard Maxwell evolution equations are satisfied. Thus 
numerical error can be expected to lead to exponential growth of the constraint for a 
hyperboloidal foliation with K > 0. The situation is more complicated in the nonlinear 
gravitational case but similar instabilities of the system of equations governing the 
constraints are associated with the extrinsic curvature |30| . A negative value of K 
(the expanding case) tends to damp constraint violation whereas a positive value (the 
collapsing case) can trigger constraint violating instabilities. 

In the collapsing direction, we perform the runs with a harmonic time slicing 
to prolong the approach to the singularity, as previously done by Garfinkle |32j . 
Results for the Hamiltonian constraint for the various codes are shown in Fig. [15] 
for the collapsing case. All the codes now show some growth in the Hamiltonian 
constraint, either of a slow or runaway type. The slow growth, exhibited for example by 
the AbigeLharm, AELCactusEinsteinADM and KrancNOR codes, can be attributed 
to the analytic constraint instabilities discussed in [30j : the Hamiltonian constraint 
violation remains small (« 10~^) at the end of the run. The runaway growth exhibited 
by the LazEv_BSSN code can again be controlled by numerical dissipation. This 
is demonstrated by the convergence results shown in Fig. [TH] for the fourth order 
dissipated version of the code. 
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Figure 15. Comparison plot of the Loo norm of tiie Hamiltonian constraint vs 
liarmonic time t for the polarized Gowdy Wave evolved in the collapsing direction 
with the p = 4 resolution. 




Figure 16. Convergence results for the Loo norm of the Hamiltonian constraint 
for the polarized Gowdy Wave evolved in the collapsing direction by the 
4"'order system LazEv_BSSN code with dissipation. After rescaling the results 
for the p = 2 by 1/16, they closely match those for the p = 4 resolution. The 
figure shows stability and convergence of the Hamiltonian constraint up to 1000 
crossing times and demonstrates good performance of the LazEv_BSSN code if 
dissipation is added. 
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The choice of specifications given in |Appcndix A.5| provides a Gowdy testbed 
capable of good discrimination between different formulations. Results for both the 
expanding (Fig. [T^ and collapsing (Fig. [T5)) directions show a wide spread in the 
performance of the different codes. We observe, as in the gauge wave test, that the 
BSSN-based codes have less satisfactory performance. 

8. Conclusions 

This first round of tests, although modest in scope is a good start at establishing the 
methods for code verification that have been deemed necessary for any complicated 
computational discipline, such as numerical relativity, to fulfill its scientific potential. 
As observed by Post and Votta [3^ in their study of the verification and validification 
of large scale computational projects, "the peer review process in computational 
science generally doesn't provide as effective a filter as it does for experiment or theory. 
Many things that a referee cannot detect could be wrong with a computational science 
paper. . . The few existing studies of error levels in scientific computer codes indicate 
that the defect rate is about seven faults per 1000 lines of Fortran" . Their observations 
are especially pertinent for numerical relativity where validation by agreement with 
experiment is not available. 

Several problems have been encountered in the course of this project. One 
problem was getting prompt response from a broad set of groups with many other 
pressures. The Apples with Apples workshops were very successful in this regard 
and were absolutely essential in jump-starting and continuing the project. But after 
the participants dispersed from the workshops, outside pressures led to predictable 
difficulties. Besides teaching and administrative duties, the overriding scientific 
pressure in the field has been solving the two black hole problem and supplying 
waveforms. This raises a complicated juggling of priorities between black hole 
simulations and code verification. In order for code verification to be attractive, the 
tests have to be useful and the investment in time has to be minimal. This adds 
emphasis on the need for tests that are simple to carry out and simple to document 
the results. 

Another level of complication in this project arises from the feedback between test 
design and the analysis of test output. This has led us to improvements in the tests and 
to their better understanding. In the robust stability test the correct interpretation of 
results for weakly hyperbolic algorithms required rethinking the proper choice of norm 
and refinement procedure for judging stability. In the gauge wave tests, the desire for 
computational efficiency in detecting nonlinear problems at an early time has led us 
to the adoption of a higher amplitude ^ = 0.5 for the test, as opposed to the original 
specifications A = 0.01 and A = 0.1. 

The robust stability test is presented as a pass/fail test. For the linear wave 
test the amplitude and phase errors in the output data for the wave profile provide 
a good comparison of code performance. For the gauge wave and shifted gauge wave 
tests, a prime challenge is the suppression of long wavelength nonlinear instabilities in 
the analytic problem. For the Gowdy test, there were unanticipated shortcomings in 
the output content that should lend valuable experience in the design of future black 
holes tests. Useful benchmarks have been established for the linear wave, gauge wave, 
and Gowdy wave tests, which have revealed clear deficiencies in various codes. Such 
deficiencies raise a clear alert that it is necessary to apply or rechcck other verification 
techniques, such as convergence tests. 
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These first round results provide a good basis for proposing new tests. Already, 
they have prompted addition of the shifted version of the gauge wave test, in which 
a non-vanishing shift fills a gap in the four original tests for periodic boundary 
conditions. This test has been useful in developing analytic and numerical techniques 
for controlling instabilities [23l [25]. A second round of boundary tests based 
upon the periodic tests have been proposed. The specifications are given on the 
Alliance website [5]. Results of some of these boundary tests have been reported 
elsewhere [241111]. The next stage is to formulate tests involving black holes. 

The code comparisons have proved useful for designing code improvements and 
for stimulating the use of new numerical techniques. During the course of this work, 
results of the shifted gauge wave test were key to recognizing the importance of discrete 
energy and flux conservation for harmonic code performance |23j . The need to carry 
out the tests with a wide range of formulations has led to the development of symbolic 
code generation [35]. Although the tests were designed for finite difference codes, they 
have been adapted and applied to pseudo-spectral codes [53]. Further independent 
studies based upon the tests have played a major part in thesis research [36l |37] . 

Establishment of the CVS data repository has been an important step in the 
documentation of test results. Instructions for accessing the data are given at [5]. 
The CVS directory structure has been significantly streamlined and documented since 
the beginning of the project. However, the difficulties in completing this analysis 
of the first round of tests has emphasized the need of a uniform standard for data 
structures and output. Rather than trying to anticipate a complete list of useful output 
quantities, it seems more desirable to output the 3-metric and extrinsic curvature at 
specified times. Then other output quantities can be constructed in post processing. 
Ideally, this should be done in some standardized way using automated routines and 
graphical interfaces. All of this would require considerable infrastructure to provide 
hardware for data storage and software for processing. This is one of the important 
matters that will be presented for discussion at future Alliance meetings. 
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Appendix A. Revised testbed specifications 

We present here the updated specifications for the five standardized testbeds. For 
each test we provide the 4-nietric of the spacetime, except for the robust stability test 
where only the initial Cauchy data is specified. This determines the 3-metric h^i, = 
dfiiy + nfj^rii, (where is the future directed unit normal to the Cauchy hypersurface) 
and the extrinsic curvature K^j,^. We use the convention A'^^ = —hPyi^Up for which 
the trace K is negative for an expanding cosmology. In all cases, the evolution takes 
place in a fixed rectilinear coordinate domain with periodic boundary conditions, i.e. 
a 3-torus. The identified "boundaries" in the 3-torus picture are located a half step 
from the first and last grid points along each axis. 

Even though we are concerned with 3-dimensional codes, for tests with only one- 
dimensional features in the x-dircction it is efficient to use the minimum number of 
grid points in the trivial y and z directions, i.e. to run the test in a long channel 
rather than a cube. For standard second order finite differencing this implies that we 
use 3 or 4 points in those directions. For all such ID tests, the evolution domain is 

X e [-0.5,4-0.5], y = 0, z = 0, (A.l) 

with grid 

x = -0.5+ (n- ?i = 1 . . . 50p, = 1/(50/9), p G Z.(A.2) 

(In the Gowdy wave test, the grid is aligned with the z-direction.) The 2D tests have 
evolution domain 

X e [-0.5; +0.5] y e [-0.5; +0.5], z = (A.3) 

with both the x and y grids satisfying (jA.2p . The parameter p allows for grid 
refinement. The coarsest p = 1 grid is useful only for debugging. Convergence tests 
should be made with p = 2 and p = 4, with benchmarks for norms, constraints, etc. 
provided by p = 4. 

We have dropped the original requirement that the tests be run with an iterative- 
Crank-Nicholson algorithm since Runge-Kutta time integrators have since proved to 
be more effective and have been commonly adopted. For each test, the size of the 
timestep dt is given in terms of the grid size to lie within the CFL limit for an explicit 
evolution algorithm. (For some codes this may be inappropriate and some equivalent 
choice of time step should be made.). A final time T, and intermediate times for 
data output, are specified for each test. They are chosen to incorporate all useful 
features of the test without prohibitive computational expense. Except for the robust 
stability test, it is important to calculate the convergence rate of the numerical error. 
Additional output variables might be essential to assess the performance of a particular 
formulation. 

Appendix A.l. Robust stability testbed 

The 3-metric is initialized as hij = (5ij+ey, where are independent random numbers 
at each grid point. All other evolution variables are initialized in the same way. The 
amplitude of the random noise is scaled with the grid as 

e e (-10~i7p2, +10-i7p2). (A.4) 

The range of the random numbers ensures that effects are below round-off accuracy 
so that the evolution remains in the linear domain unless instabilities arise. 
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The timcstcp is specified to be dt = dx/10 = 0.002//9. The use of 4 distinct 
gridpoints in the y and z directions aUows for instabihties associated with the 
checkerboard mode. 

The test should be run until one is confident that dissipation effects do not cloud 
the result. Without artificial dissipation, a run time of one crossing time, using output 
at every time step, is usually sufficient. This corresponds to 500p time steps. The 
test is passed if the norm satisfies the inequality (fTTjl for all resolutions, for some 
fixed choice of constants A and K. Appropriate norms for both first and second order 
systems are recommended in [TT] and are publicly available as Cactus thorns [55] . 



Appendix A. 2. Linear wave testbed 

The initial 3-metric and extrinsic curvature Kij are given by a transverse, trace-free 
perturbation with components 

ds^ = -dt^ + dx^ + {1 + H) dy^ + {1-H) dz^, (A.5) 

where 

H = Asu.{^-^^^y (A.6) 

This describes a linearized plane wave traveling in the x-direction. The wavelength is 
set to d = 1 to match the periodicity of the evolution domain. The metric has lapse 
a = 1 and shift = 0. The nontrivial components of extrinsic curvature are 

In order to test 2-dimensional effects, the rotation 

x = ^{x'~y'), y^l={x' + y'). (A.8 

leads to a wave propagating along a diagonal. The resulting metric is a function of 
'2ti{x' -y' - ty/2) 



Kyy = -7;dtH, K,,^-dtH. (A.7) 



sm 



, where d' = dV2 . (A.9) 



To obtain the required periodicity of the evolution domain, we set d = 1 in the ID 
simulation and d' = 1 in the diagonal simulation. The test should be run in both 
axis-aligned and diagonal form. 

The test is performed with amplitude A = 10~^, so that quadratic terms are of 
the order of numerical round-off. The time step is set to dt = dx/ 4 = 0.005/p As 
in the gauge wave case, the ID evolution is carried out for T = 1000 crossing times, 
i.e. 2 X 10^/5 time steps , with output every 10 crossing times. The 2D diagonal runs 
are carried out for T = 100, with output every crossing time. The output quantities 
are the Loo and L2 norms, the maxima and minima, and profiles along the x-axis 
through the center of the grid of gyy, g^z, Hamiltonian constraint; and the Loo error 
norm for g^z (measuring the difference from the exact solution). 



Implementation of standard testbeds for numerical relativity 



26 



Appendix A. 3. Gauge wave testhed 

The test is based upon the 4-mctric 

ds^ = (1 - H){~dt^ + dx^) + dy^ + dz^, (A.IO) 

with H given by (lA.Op . which is obtained from the Minkowski metric ds^ — —dt^ + 
dSp' + dif' + dz^ by the transformation 



x+ 4^ cos 

47r 



[r^] ' (A.ii) 



y ^ y, 

z — z. 



This describes a sinusoidal gauge wave of amphtude A propagating along the x-axis. 
The extrinsic curvature is 

K - ^"H^j .^12) 

Kij — otherwise. (A. 13) 

As for the linear wave, the rotation (|A.8|) leads to wave propagation along a diagonal 
with periodic boundary conditions. 

The gauge wave is run with amplitude A = .5. The time coordinate t in the 
metric (jA.lOp is harmonic and the gauge speed is the speed of light. The time step 
is set to dt = dx/A = 0.005//9. The ID evolution is carried out for T = 1000 crossing 
times, i.e. 2 x 10^ p time steps (or until the code crashes), with output every 10 crossing 
times. The 2D diagonal runs are carried out for T — 100, with output every crossing 
time. 

Output consists of the Loo and L2 norms, the maxima and minima, and profiles 
along the .T-axis through the center of the grid (y = z = 0) of gxx, ct, tr{K) and the 
Hamiltonian constraint; and the L2 error-norm for g^x- 



Appendix A. 4- The shifted gauge wave test 

The shifted gauge wave is obtained from the Minkowski metric ds^ = —dt^ + dx^ + 
dy^ + dz^ by the harmonic coordinate transformation 




(A.14) 



y = 2/, 



z — z 

which leads to the Kerr-Schild metric 

ds^ = -dt^ + dx'^ + dy'^ + dz'^ + Hk^kpdx'^dx'^ (A.15) 

where 



(A.16) 
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and H is again given by (|A.6|) . The extrinsic curvature is 

K,j = otherwise. (A. 18) 

This metric describes a shifted gauge wave of amplitude A propagating along the x- 
axis. The coordinate transformation (jA.SP rotates the propagation direction to the 
diagonal. 

The shifted gauge wave test is run in a harmonic gauge with amplitude A = 0.5 
in both ID form and diagonal 2D form. As in the linear wave test, for the required 
periodicity we set d = 1 in the ID simulations and d' = 1 in the 2D simulations. We 
set the timestep dt = dx/4 = 0.005/p. The ID evolution is carried out for T = 1000 
crossing times, i.e. 2 x 10^ p time steps (or until the code crash). The 2D runs are 
carried out for T = 100. 

Output data consist of the profiles along the .x-axis through the center of the 
grid (y = z = 0) of gtt, gxt-, and gxx, the L2 and L^o norms of the error and of the 
Hamiltonian constraint. 

Appendix A. 5. Polarized Gowdy wave testbed 

The polarized Gowdy metrics describe an expanding, toroidal universe containing 
plane polarized gravitational waves with metric 

ds^ = t-^/^e^/\-dt^ + dz^) + te^dx^ + e-^dy\ (A.19) 

where A and P are functions of z and t only and are periodic in z. The universe 
expands as t increases. The test is carried out in both the collapsing and expanding 
situations. The metric is singular at t = 0. 

The Einstein equations reduce to a single evolution equation 

Pjt+t-^P^t-Pzz^O (A.20) 

and the constraint equations 

A,t = t(Pj + p2) (A.21) 

and 

A,, = 2tP,Pt. (A.22) 
The test is based upon the particular solution to (jA.20[) 

P = Jo{2TTt) cos(27rz), (A.23) 
where J„ are Bessel functions. The metric and extrinsic curvature are 

9xx = te^, gyy - te-P, g,, = t-^'^e^/'\ (A.24) 



h^/^e-^/^e^il+tP,)^ 



2 



Kyy-~\t^'^e-^'^e--il-tP,l (A.25) 



Vi/v/4(t-i-A,o, 



with 



trK = -^t^/*e-^^\3t'^ + A,t). (A.26) 
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The shift vanishes and the lapse is 

a = = t'^/'e'/\ (A.27) 

For the choice (iX23|) . the constraints (|A.21|A.22p yield 

A = -27rtJo(27rt)Ji(27rt)cos2(27rz) + 27r2i2[j2(27rf) + J2(27rt)] , 
- ^{i27r)^[j^{2n) + JH27T)] - 2^Jo(27r) Ji(2^)}. 

While P slowly decays to zero, A undergoes linear growth due to the cosmological 
expansion, and both P and A exhibit gravitational wave oscillations. 

The velocity of light is constant in the coordinates chosen in (jA.lQp so that, 
with a fixed spatial discretization dz, the Courant condition is consistent with a fixed 
timestep dt. This makes the gauge (|A.f9P convenient for evolving in the expanding 
direction by choosing the initial data from the exact solution at t = 1, which yields 
data of order unity. 

In the backward in time evolution, wc choose a harmonic time slicing r which only 
asymptotically reaches the singularity. Starting with the metric (|A.19|) . the slicing is 
obtained by a transformation t — F(t), where the harmonic condition □ r = implies 
F{t) ~ kef^ . In order to start the collapse slowly, the free constants c and k arc chosen 
so that the new lapse satisfies a = 1 at the initial time t = Iq. This is accomplished by 
picking to for which Jo(27rto) = so that (|A.28p implies a is independent of z. Using 

To = -lnf^V A(fce^^°,z)-Ao 



c \k , 
we obtain 

ao = ctl^^ e^°^\ (A.29) 

Given our requirement cio = 1, and choosing to ~ to, i.e. F(to) = tq, we get 

c = tg^/* e-^«/^ k = toe-'''". (A.30) 

Wc choose a particular value of to such that the initial slice is far from the cosmological 
singularity, but not so far that wc have to deal with extremely large numbers by picking 
the 20th zero of the Bessel function Jo(27rto), which yields to ~ 9.8753205829098, 
corresponding to 

c ~ 0.0021195119214617, k - 9.6707698127638. 

The time step is set to dt ^ dz / 4 ^ 0.005/p with run time T — 1000 or until code 
crash. Output consists of the Loo and ij2-norms, the maxima and minima, and profiles 
along the z-axis through the center of the grid of gzz, ct, tr{K) and the Hamiltonian 
constraint. We output norms every crossing time, and profiles either every 10 crossing 
times or once per crossing time, depending on the behavior of the simulation. We also 
output the Loo error norms of the difference from the exact solution for g^x and gzz 
for the expanding direction. 



Appendix B. Code descriptions 

Appendix B.l. Standard ADM: Kranc_FreeADM, and AELCactusEinsteinADM 
codes 

The formulation of the Einstein equation by Arnowitt, Deser and Misner (ADM) [33] 
provides a standard notion for "evolving" space-time as an initial value problem in 
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general relativity, which was initially presented in a Hamiltonian context. What is 
referred to as a "standard ADM" system in the numerical relativity community is 
a reformulation due to York [JHIj which one obtains by 3+1-decomposition of the 
Einstein tensor (as opposed to 3+1-decomposition of the Ricci tensor in the original 
ADM version) , or equivalently by adding appropriate constraint terms to the evolution 
equations. As pointed out by Frittelli [41], York's "standard ADM" system does in 
particular have nicer properties regarding the constraint propagation system. This 
system is particularly simple, has a long history in numerical relativity and exhibits 
some typical problems. We therefore use it as the starting point for our numerical 
comparisons. The evolution equations are 

daij = ~ 2aIUj + V,/3j + VjA (B.l) 

dtK,, = aRf + aKK.j - 2aK,kK''j - V,V,a 

+ (V,/3'=)/u, + (Vj/3'=)Xfe, + (3^VkK,j, (B.2) 

and the constraint equations are 

j^^^ADM ;= Ri^) + - KijK'\ (B.3) 

Mi = Mf^"' A'^, - V,A:, (B.4) 

where {"jij^Kij) are the induced three-metric and the extrinsic curvature, {a,f3i) are 
the lapse function and the shift covector, V, is the 3-dimensional covariant derivative 

(3) 

and R^j is the 3-dimensional Ricci tensor associated with jij . 

We have tested two implementations of the standard ADM system, the 
code AELCactusEinsteinADM, which is freely available via the website |42], and 
Kranc-FreeADM which is based on the Cactus Toolkit [42] and Kranc software 
[35] . AELCactusEinsteinADM uses a hardcoded ICN time update scheme (see e.g. 
[TT|). whereas Kranc_FrceADM uses a method of lines (MoL) approach based on the 
CactusMoL thorn (in practice, RK3, RK4 and ICN (see e.g. [11]) have also been 
used, as indicated). In all of these codes, spatial partial derivatives are reduced to 
partial derivatives of the 3-metric, i.e., all expressions such as Christoffel symbols 
are expanded out. Due to the absence of first-order variables, no further ambiguities 
arise. Centered second and fourth order discretization is used (see [Appendix C3] ) , 
and third order Kreiss-Oliger dissipation is optionally applied to all variables (see 
[Appendix C.2[ ). 

The hyperbolicity of the ADM free evolution scheme has been analyzed and found 
to be weakly hyperbolic with the type of gauge conditions that we use |11] . Since many 
of our tests are essentially ID tests, where ADM yields good results, we have also 
analyzed the hyperbolicity of ADM in ID. For simplicity of presentation we restrict 
ourselves to the linearized case. Assuming propagation in the x-direction we obtain 
the following evolution equations. For the off-diagonal components. 

The evolution equations for and ^xz are analogous to the evolution equation for 
7y2. The fact that the evolution equations for K^y and K^z are trivial renders the 
evolution system for the off-diagonal components weakly hyperbolic, see e.g. [TT]. For 
the diagonal components, 

dtlii ^'i.Kii {i = x,y,z), (B.5) 

dtKxx = dxxOi + \dxx{lyy + Jzz), (B.6) 
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dtKjj = -d^^-fjj {j = y,z). (B.7) 

Considering for simplicity the densitized lapse case, a = ^/j, the evolution 
equation for K^^ becomes 

dtKa;x = ]^dxxlxx + dxxhyy + Izz) 

and one finds that the diagonal subsystem is only weakly hyperbolic. However, within 
the subclasses of gauge wave (7^^ = = 0) or linear wave {'^xx = 0) data, the ID 
ADM system corresponds to copies of the ID wave equation and is therefore well- 
posed. 



Appendix B.2. AhigeLharm 



The Abigel code developed in Pittsburgh is based upon a symmetric hyperbolic 
formulation of the Einstein equations using generalized harmonic coordinates 
satisfying the curved space wave equation 



1 



1 



(B., 



\/-9 V~5 

where H°' are harmonic source terms. The original version of the evolution equations 
was [35] 

r'^do.dpr'' = S^" (B.9) 

where the left hand side is the principle part and the right hand side contains 
nonhnear first-derivative terms. Here g'^" = sj—gg^", with g = <iet{g^u) = det(5'"'). 
and the harmonic constraints dug'^" = iJ^ are used in the Einstein equations to 
eliminate second derivatives in the source terms 5'"^. For further details concerning 
the formulation and its implementation see [43j . 

The code with which the tests were performed was constructed by rewriting (jB.9[) 
in the flux conservative form 

da {g'^^dpr") = S^'"- 
and reducing it to the first order in time form 



9 



-9, 



9 



9- 



9*'9'' 



5" 



in terms of the evolution variables (5^'', Q'^"), where 



B.IO) 

B.ll) 
B.12) 
B.13) 
B.14) 

B.15) 
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and /i'-' = g*-' — g''^g-'^ /g^^ is the spatial 3- metric. Centered derivatives are used to 
finite difference (jB.lip and the source terms S^'^ in (|B.14p . The remaining part of 
Eq. (|B.14p is finite-differenced as foUows: 

9[I+1/2,J,K] - I — 2 j 

r,ti tj 



h 



ij _ ij _ 9[I+1/2,J,K] 9[i+i/2,J,K] 

[1+1/2, J,K] ^ 9[I+i/2,J,K] „tt 



9[I+1/2,J,K] 



B.16) 

B.17) 
B.18) 
B.19) 

\9 ''[I,J,K] \9[I+l/2..LK] / 

where the averaging operator A+x is defined in [Appendix C.l[ The code is evolved 
as a first differential order in time and second order in space system with a 2-step 
iterated Crank-Nicholson algorithm or 4th order Runge-Kutta integrator. 



[h-^dyr"^ ^^^-D^x (/^f/+i/2,xKi ^+./?o.CAi^]) + '^(^') 



9''* r^^A _n f 9lI+l/2,.LK] ^ 



Appendix B.3. HarmNaive 

The HarmNaive code is based upon harmonic coordinates but differs from the 
AbigeLharm code because the evolution system consists of only the 6 wave equations 
(jB.lOp for the spatial components g*-'. The time components are propagated by the 
harmonic conditions (jB.8[) . i.e. 

atr* + a,r' = ^"- (B.21) 

The coupling between g^^ and g"* makes the system only weakly hyperbolic. 

The evolution equations for g^^ and Q'-' are finite differenced as in the AbigelJiarm 
code. The evolution equation (|B.2ip for g"* is approximated by central differences. 
The update scheme is a 2-step iterative Crank-Nicholson algorithm. 



Appendix B.4- KrancNOR code 

Appendix B.^.l. Continuum formulation: Nagy, Ortiz and Reula suggested [17] 
modifications to the ADM system such that it can be made strongly hyperbolic whilst 
remaining in second order form. The system we use includes slight adjustments of [9]. 
Additionally, we use an evolved lapse. 
The variable /; is defined as 

f^^l''\l^k^-\plkLi) (B.22) 
with parameter p. This introduces the new constraint Gi where 

G^ /. - l''\l^k,l - \piku)- (B.23) 

Starting from the ADM evolution equations, an evolution equation for /; is 
obtained by differentiating (|B.22p and commuting space and time derivatives. The 
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Hamiltonian and momentum constraints are added with parameters c and b, and 
derivatives of the Gi are added with parameters a and a': 

dtlij = - 2aKij 

dtK^j = - D^Dja + - 2A',fci^^- + K^jK) + ^G(,,,) + {cH + a'Gk^n'^'H, 

= -aF{a,K,x'). 

The variables 7^^, iiTy , fi and a are evolved. Due to the symmetries of 7y- and Kij, 
this leads to 16 evolved variables. We write the Ricci tensor entirely in terms of 7^-; 
fi is only used where it appears as part of G;. 

For those tests requiring harmonic slicing, the lapse source function is 

F{a,K,x')=aK (B.24) 

and for the expanding Gowdy test, 

Fia,K,x')=K33/a (B.25) 

which is compatible with the exact lapse in this case. We make the following choice 
of parameters: 

a=l, 6^=1, a' = 0, p==2/3, c = 0. (B.26) 
Note that choosing parameters 

a = 0, 5 = 0, a' = 0, p = 0, c (B.27) 
leads to a standard ADM system. This is useful for testing the code. 

Appendix B.^.2. Semi-discrete scheme: To form the semi-discrete approximation, 
discretization in space is performed according to the standard second order accurate 
discretization IC.ll 

Finite differences are taken only of the evolved variables 7^^-, Kij, fi and a. This 
means that where derivatives of other quantities appear, they are explicitly written in 
terms of derivatives of the evolved variables (e.g. by using the Leibniz rule). 

We do not add Kreiss-Oliger type artificial dissipation, as it was not necessary 
for stability. 

Appendix B.^.3. Time integration: Time integration is performed using the method 
of lines with the iterative Crank-Nicholson (ICN) method. 

Appendix B.4-4- Output: For our state vector v = {jij , Kij , fi)'^ we define the L2 
and norms: 

= iv'^v'^vlki + v'^v'^K.^Kki + v''nfj)h' (B.28) 
grid 

Ml+ ^ Ml, + E (^*S^'S""i^+™7.yi^+n7fci)/^' (B.29) 
grid 

where ry = diag(l, 1, 1). This is the norm obtained from a reduction to first order 

of the semi-discrete equations. The exact solution is denoted u" = u{t",Xj) and the 

error is defined as 

£ = v-u. (B.30) 
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For the stability test, the exact solution is taken to be Minkowski in Cartesian 
coordinates. For those tests which are perturbations of this solution, we analyze 
relative error with respect to this background. We denote the background solution as 
Ub- Hence the relative error about this background is 

_ II^IU. 



\\u - ub\\l2 

In general, we run mitil this quantity exceeds 0.2 (a relative error of 20%). 



(B.31) 



Appendix B.5. Family of BSSN (Shihata-Nakamura and Baumgarte- Shapiro) 
formulations 

The family of BSSN systems is constituted by variations of an evolution system 
that had originally been proposed by Nakamura in the late 80s, and has been 
subsequently modified by Nakamura-Oohara and Shibata- Nakamura [HlllSlllS], and 
later by various other authors. The formulation is characterized by introducing a 
contracted connection term as a new variable, a conformal decomposition of the 
metric and extrinsic curvature variables, and adding constraints to the evolution 
equations. In particular, the system can be viewed as the NOR-system plus a 
conformal decomposition which leads to the evolution of a unimodular metric. The 
advantage of this formulation was re-announced by Baumgarte and Shapiro |47| . 

Modifications of the system have been obtained by variations in how derivatives of 
the new variables are written, how the gauge is specified, how algebraic constraints are 
treated, and the way (differential or algebraic) constraints are added to the evolution 
equations. A detailed discussion of well-posedness for the BSSN family has been given 
by Gundlach and Martin-Garcia [51 [HI [IHIj to which we refer for details about the 
BSSN family. 

The set of evolved variables are the logarithm of the conformal factor the 
conformally rescaled three- metric 7ij, the trace of the extrinsic curvature if, the 
conformally rescaled traceless extrinsic curvature Aij, and the contracted Christoffel 
symbols F*: 



<^ = (l/12)log(det7„), (B.32) 

1^J -e-'^^7,„ (B.33) 

K ^Y'K,,, (B.34) 

A,j = e'^^{K,j - (l/3)7,,if), (B.35) 

r =fj,7^^ (B.36) 

This immediately leads to the two algebraic constraints 

det7y=l, i* = (B.37) 
and the differential constraint 

f^ _ ^jfcf ^ 0, (B.38) 



which are9 propagated by the evolution equations. Note that densitized quantities 
(those with a tilde) have their indices raised and lowered with the conformally rescaled 
three- metric 7^^ . 
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The standard Hamiltonian and momentum constraints (|B.3IB.4p and (|B.4p can 
be expressed in the BSSN variables as 

- iyi'^ - {2/2,)AK, (B.39) 

M^ = QAU{bj^) - 2^(Av) - {2/i){D^K) + 7^^(^jifc,;). (B.40) 

The BSSN evolution equations, which are obtained from the ADM equations (jB.ll 
- IB.4P by using the definitions (|B.32I - lB.36p and making a standard choice for adding 



constraints, are 

£„(^ = - {l/&)aK, (B.41) 

Cnl^J = -2aA,j, (B.42) 

Cr,K = ~ D'D.a + aA.jA'^ + {l/?,)aK^, (B.43) 

£„1,^. = - e-^^^iD^Djaf^ + e-4'^a(i?,g^'5^)^^ + aKA,, - 2al,fci^■, (B.44) 

£„r = -2(9j-a)i*^'+2a(f}^^i'^^-(2/3)fJ(ajA') + 6i*^'(aj(p)), (B.45) 
where Di is the covariant derivative associated with 7^^, and £„ ~ dt ^ Cp is the Lie 



derivative along the miit normal. Note that J CnKd^x is positive definite apart from 
boundary terms involving the lapse (which vanish for periodic boundary conditions). 
The Ricci curvature Rf/'^^ in terms of the BSSN variables becomes 

nBSSN _ fy , nV> 

J^ij -fly -f ^ijj 

R,, = - {l/2W'^l^kl^J + 7fe(.5,)f^ + f^f + 27'"f f(,f,)fc^ + 7'"f Lf fei,. 

Note that there are different ways to numerically compute the trace free part of the 
Ricci tensor, e.g. one can project out the trace of the Ricci tensor according to 

Rf/ = R,, - ii?7.,, (B.46) 

compute the Ricci Scalar from the Hamiltonian constraint (|B.39p . or compute the 
trace free part explicitly by assuming the algebraic constraints hold. 

We refer to the code descriptions below for details concerning the individual codes. 
In summary, the fundamental dynamical variables in BSSN are ((^,7^, KAijA^)^ 
which total 17. The 4 gauge quantities are (a,/?*). 

Appendix B.5A. Concrete implementations We have compared a number of codes 
based on variants of the BSSN system. Several of these are based on the Cactus 
computational toolkit [42]: the CCATIE.BSSN Hi |49] and Kranc_BSSN [50] codes, 
and the LazEv.BSSN [SJ code. Of these, CCATIE_BSSN and Kranc_BSSN use the 
CactusMoL time integrator, which provides the RK3, RK4 and ICN methods, among 
others (see e.g. [H]). KrancJBSSN is based on the Kranc code generation software 
package [55] . 

All codes use straightforward replacement of partial derivatives by standard 
second order centered finite differences with a three point stencil (most codes are 
also able to use standard centered fourth order finite differencing). 

Most of the BSSN codes have a long history of use in production environments 
and have a large number of parameters that allow them great fiexibility, e.g. 
regarding details of the numerical methods, gauge conditions, or the way the algebraic 
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constraints are treated. Typical options to solve the algebraic constraints at every 
intermediate timcstcp use the following replacements: 

• Ensure that 7,., has unit determinant by setting 

• Ensure that Aij remains trace-free by setting 

1 ~ 
3 

Divide Aij by the same factor that is used to remove the determinant of 7^ : 



A, - -in - jAm7''7'"- (B.48) 



Note that an ambiguity arises whenever F* or j^'fij.k occur, as they are related 
analytically by the equation F* = ,i — ^7''(ln7),;. If the constraint 7 = 1 holds, 
e.g. if it is enforced at each timestep, this is equivalent numerically (up to round-off 
error) to F* = — 7'-^ ,j . Some authors replace 7*-' ,j using — F* only when the expression 
appears under a derivative, but more complicated rules have also been applied. 

Ref. |52] describes a widely used combination of BSSN system and gauge condition 
in detail and examines this system's hyperbolicity. 



Appendix B.6. KrancFN 

Appendix B.6.1. Continuum formulation: The Friedrich-Nagy system |53j is a frame- 
based first order formulation that has been shown to yield a well-posed initial boundary 
value problem. The formulation starts from the four dimensional vacuum equations 

:=[e,,e,/]''-(F7^j-Fj^/)e;f'' = 0, ^ = 0,1,2,3 (B.50) 

A/jKL := Rijkl{T) - CijKL = (B.51) 

HjKL := V/C/xl' = 0, 7 = 0,1,2,3 (B.52) 

where e/ denote the tetrad vectors with coordinate components e/^; and Tj^ j are 
the connection coefficients defined by Ve^eK = F/'^ifCj and satisfying ijjm ^ i'^ k + 
rjKj^i'^M = 0. RijKL and Cijkl denote the components of the Riemann and 
Weyl tensor with respect to the tetrad. The Riemann tensor is given in terms of the 
connection coefficients by 

Rij'-KiT) = e/(Fj^K) - ej(F/K) 

M I J ^ i- I J M + ^ M J / + t / J K- (13.^6) 

Equation (|B.50P states that the connection is torsion free, ()B.5ip are the vacuum 
Einstein equations and (jB.52p is the Bianchi identity for a vacuum spacetime. From 
(|B.50[) - (jB.52[) . a symmetric hyperbolic evolution system is obtained by choosing 
certain combinations of components of the above equations as well as a gauge that is 
adapted to the boundary. 

Assuming a boundary at z = const, we foliate the interior domain by time-like 
hypersurfaces Tc given by z = c = const. The frame is adapted to this foliation and 
boundary such that the frame vector 63 is orthogonal to Tc, which implies for the 
coordinate components 

Ca^^O, a = 0,1, 2, > 0. (B.54) 
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63 being the unit normal to implies F^'^h = r^^^^). 

The mean extrinsic curvature of T^, is prescribed as a function of the coordinates 
f{x^) and used to eliminate the connection coefficient Fq'^o from the equations, 

T^\ = f + Ti\+T2\. (B.55) 

The variation of Cq within is prescribed by functions F^[x^)^ A = 1,2 according 
to DggCo = F'^eA, where D denotes the induced connection on T^. This eliminates 
the connection coefficients 

Tq^o=F'^, ^=1,2. (B.56) 

The tetrad vectors arc Fermi-transported along cq with respect to D and therefore 

Fo^B = 0, A,B^ 1,2. (B.57) 

The coordinates {x^} are chosen such that the tetrad vector cq represents the time 
flow dt, i.e., 

eo^' = (B.58) 

The ten independent components of the Weyl tensor are encoded in the symmetric 
and tracefree tensor fields 

Eij := Ciojo, Bij := —Coiki^''''£''^ j 

corresponding to the electric and magnetic parts with respect to cq. The conditions 
S^^ Eij = 5^^ Bij = are incorporated explicitly by eliminating 

E33 = -(-Ell + E22), B33 = -{Bn + B22) (B.59) 

from the equations. In total the Friedrich-Nagy system has 37 variables, namely 

u= (eA^e3^,F,o„F3%,F(^3^),FA^c,£;.A,B»A)^, (B.60) 

where 

A,B,C^1,2, i,j = l,2,3, p = 0,l,2, ^ = 0,1,2,3. 

A symmetric hyperbolic evolution system for the variables (|B.60[) is obtained by 
taking the following combinations of (|B.50p - (jB.52p : 

W = 0, ro3^ = 0, AoBab^O, Aoi31 = 0, Ao232 = 0, 

A0132 + A0231 = 0, A0130 + A1232 ~ 0, A0230 + A2131 — 0, 
A^B03 = 0, A^ooa = 0, A3A03 + A303A — 0, ?7"''A3at3 ~ 0, 

where the convention for the indices is the same as in Eq. (jB.60[) and a, 6 = 0, 1,2. 
The resulting system is given explicitly in [53l [36] and is of the form 

A°atu + A*(9,u + B(u,F) = 0, (B.61) 

where F = [f , F^ ,d^f ,d^F^) represents the gauge source functions and their 
derivatives. The matrices A°,A* are symmetric and depend on the coordinate 
components of the frame. A" is positive definite as long as 1 — (ei°)^ — (e2°)^ — (63°)^ > 
0, which corresponds to eo being time-like. Characteristics are time-like and null. 
The remaining components of (jB.50|) - (jB.52[) . 

T,/ = 0, A,/k-0, ii/jfcoe^'Vn = 0, 



TJ A3 ^ 3' W ^mn r. ^ tt ,mk , , 3m rr n 
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only contain derivatives in directions orthogonal to eg and are satisfied if satisfied 
initially by virtue of the evolution equations (see [Mj). Cq in general is not hypersurface 
orthogonal and therefore the constraints do contain derivatives in direction of dt- In 
order to monitor these constraints during a numerical evolution, we eliminate the time 
derivatives by means of the evolution equations. 



Appendix B.6.2. Numerical implementation: The code is based on the Cactus 
Computational Toolkit [H] and the Kranc software [351 ISS] • The spatial discretization 
of (|B.61[) is done in a straight forward way 

dtn = -(A°)-iA*Au + (A°)-iB(u, F), (B.62) 

where Di is the 2nd (or 4th order) accurate centered derivative operator in the 
direction i (see [Appendix G.l\ . Time integration is done with the method of lines 
(CactusMoL) using ICN for the 2nd order scheme and RK4 for the 4th order scheme. 
If needed, artificial dissipation is added to the right hand side of equation (jB.62p in 
the form 

(A")-iQrfU, (B.63) 

where Qd is the Kreiss-Oliger dissipation operator (see |Appendix C.2| . Respecting the 
symmetrizer in the dissipation term is essential; replacing it by the identity matrix 
triggered exponentially growing continuum modes e.g. for the gauge wave testbed with 
non-linear amplitude. 



Appendix B. 7. LSU-HyperGR 

This symmetric hyperbolic first order formulation is described by Sarbach and Tiglio 
in [SI]. The system has 34 evolved variables which are the standard ADM metric 
7ij, extrinsic curvature Kij and lapse a, as well as extra variables di^ij ~ dk"fij and 
Ai = dia/a, introduced to make the formulation first order in space. 

In addition to the Hamiltonian constraint Ti. and the momentum constraint A4i, 
the constraints arising from those new variables are 

Ca. ^A,-d,a/a, (B.64) 

C'kij =dkij-dkjij, (B.65) 

Cikij = d[idk]jk- (B.66) 

The system of PDEs resulting from the standard ADM 3+1 decomposition of the 
Einstein equations is only weakly hyperbolic. To get a symmetric hyperbolic system 
the principal part has to be modified further. This is done by adding the constraints to 
the right hand sides of the evolution equations with appropriate multiplicative factors 
Ci£,iViX E^nd L. Here these parameters are chosen to be constant in space, although in 
general this is not necessary. The full set of equations is then 

9o7y = -2i^y, (B.67) 

doK,^ = R,, - ^V,;V,« - 2K,aK''^ + KK,j + ^7,, W + C7"'C,(,,)b, (B.68) 

dodkij = - 2dkK,j - 2AkKij + vik{iMj) + xitjMk, (B.69) 
^oa = - F{a,K,x^') + Sixf"), (B.70) 
dF{a,K,xn , ldF{aJ<^ l dF{a,K,x^) , 
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where do = {dt — Cp)/a, Rij is the Ricci tensor and K the trace of the extrinsic 
curvature. The functions F{a, K^x^) and S{x'') are pure gauge and can be chosen 
freely. The choices 5* = and F = aK provides harmonic gauge conditions. 
Restriction of the parameters Xj ? j d to the family 

L = -1/2, Cr, - -2, e = -l/2x + 1/477 - 1/2 (B.72) 

results in a strongly hyperbolic system. A symmetric hyperbolic subfamily is given by 
= — 1, which leaves x s-s the single free parameter (constrained only by the condition 
X 7^ 0). The runs presented here were done with the specific choice of x = ~1- 

To ensure a numerically stable discretization based on the energy method for 
hyperbolic equations, second order spatial differencing operators that satisfy the 
summation by parts (SEP) condition are used [551 156] . 

Furthermore a small amount of dissipation (standard Kreiss-Oliger dissipation 
operators) is added to the right hand sides of the evolution equations. 

The integration in time is done with a third order Runge-Kutta scheme. 



Appendix C. Numerical methods 

Appendix C.l. Spatial discretization 

Most of our numerical results are based on second order accurate centered 
discretization: 

a.^D,.. s.3,^{'^^^ ^It] . (CD 



where 



DnV 



0"3 



Ax 

Ax ' 
2Ax 



D,D^., := (C.3) 

For a summary of definitions and results for standard fourth order discretizations we 
again refer to , where in particular some results concerning the evolution systems 
considered here are derived. 

Finally, averaging operators A± are defined as: 

(C.4) 

A_u, (C.5) 



Appendix C.2. Artificial Dissipation 

For second order accurate codes, it is common practice to add third order accurate 
Kreiss-Oliger dissipation [57j to all right-hand-sides of the time evolution equations 
as 



9fU — > 9fU + Qu. 



(C.6) 
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Here we use the following general form of the Kreiss-Oliger dissipation operator Q of 
order 2r, 

Q = a(-l)'-/i2'-i(£i+)V(Li_)722^ (C.7) 

for a 2r — 2 accurate scheme, where the parameter a regulates the strength of the 
dissipation and p is a weighting function, which is typically set to 1 in the interior 
but may go to at a boundary. Since we mostly focus on second order accurate codes 
here, the relevant case is r = 2, for which 

Q = -ah^{D+)^p{D_f/l6, (C.8) 

which may be implemented using Erik Schnetter's Cactus thorn AEITh.orns/Dissipation 

m- 
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